home *** CD-ROM | disk | FTP | other *** search
/ Aminet 1 (Walnut Creek) / Aminet - June 1993 [Walnut Creek].iso / usenet / sources / volume2 / aplictns / matlab / src.4 < prev    next >
Internet Message Format  |  1988-11-02  |  49KB
  1. Path: xanth!nic.MR.NET!hal!cwjcc!mailrus!ulowell!page
  2. From: page@swan.ulowell.edu (Bob Page)
  3. Newsgroups: comp.sources.amiga
  4. Subject: v02i044:  matlab - matrix laboratory, Part04/11
  5. Message-ID: <10019@swan.ulowell.edu>
  6. Date: 2 Nov 88 21:43:07 GMT
  7. Organization: University of Lowell, Computer Science Dept.
  8. Lines: 1220
  9. Approved: page@swan.ulowell.edu
  10.  
  11. Submitted-by: strovink%galaxy-43@afit-ab.arpa (Mark A. Strovink)
  12. Posting-number: Volume 2, Issue 44
  13. Archive-name: applications/matlab/src.4
  14.  
  15. #    This is a shell archive.
  16. #    Remove everything above and including the cut line.
  17. #    Then run the rest of the file through sh.
  18. #----cut here-----cut here-----cut here-----cut here----#
  19. #!/bin/sh
  20. # shar:    Shell Archiver
  21. #    Run the following text with /bin/sh to create:
  22. #    src-4
  23. # This archive created: Wed Nov  2 16:20:28 1988
  24. cat << \SHAR_EOF > src-4
  25.       IF (OP .GT. DOT) GO TO 70           
  26. C                      
  27. C     ADDITION         
  28.    01 IF (M .LT. 0) GO TO 50              
  29.       IF (M2 .LT. 0) GO TO 52             
  30.       IF (M .NE. M2) CALL ERROR(8)        
  31.       IF (ERR .GT. 0) RETURN              
  32.       IF (N .NE. N2) CALL ERROR(8)        
  33.       IF (ERR .GT. 0) RETURN              
  34.       CALL WAXPY(M*N,1.0D0,0.0D0,STKR(L2),STKI(L2),1,        
  35.      $            STKR(L),STKI(L),1)      
  36.       GO TO 99         
  37. C                      
  38. C     SUBTRACTION      
  39.    03 IF (M .LT. 0) GO TO 54              
  40.       IF (M2 .LT. 0) GO TO 56             
  41.       IF (M .NE. M2) CALL ERROR(9)        
  42.       IF (ERR .GT. 0) RETURN              
  43.       IF (N .NE. N2) CALL ERROR(9)        
  44.       IF (ERR .GT. 0) RETURN              
  45.       CALL WAXPY(M*N,-1.0D0,0.0D0,STKR(L2),STKI(L2),1,       
  46.      $            STKR(L),STKI(L),1)      
  47.       GO TO 99         
  48. C                      
  49. C     MULTIPLICATION   
  50.    05 IF (M2*M2*N2 .EQ. 1) GO TO 10       
  51.       IF (M*N .EQ. 1) GO TO 11            
  52.       IF (M2*N2 .EQ. 1) GO TO 10          
  53.       IF (N .NE. M2) CALL ERROR(10)       
  54.       IF (ERR .GT. 0) RETURN              
  55.       MN = M*N2        
  56.       LL = L + MN      
  57.       ERR = LL+M*N+M2*N2 - LSTK(BOT)      
  58.       IF (ERR .GT. 0) CALL ERROR(17)      
  59.       IF (ERR .GT. 0) RETURN              
  60.       CALL WCOPY(M*N+M2*N2,STKR(L),STKI(L),-1,STKR(LL),STKI(LL),-1)             
  61.       DO 08 J = 1, N2                     
  62.       DO 08 I = 1, M   
  63.         K1 = L + MN + (I-1)               
  64.         K2 = L2 + MN + (J-1)*M2           
  65.         K = L + (I-1) + (J-1)*M           
  66.         STKR(K) = WDOTUR(N,STKR(K1),STKI(K1),M,STKR(K2),STKI(K2),1)             
  67.         STKI(K) = WDOTUI(N,STKR(K1),STKI(K1),M,STKR(K2),STKI(K2),1)             
  68.    08 CONTINUE         
  69.       NSTK(TOP) = N2   
  70.       GO TO 99         
  71. C                      
  72. C     MULTIPLICATION BY SCALAR            
  73.    10 SR = STKR(L2)    
  74.       SI = STKI(L2)    
  75.       L1 = L           
  76.       GO TO 13         
  77.    11 SR = STKR(L)     
  78.       SI = STKI(L)     
  79.       L1 = L+1         
  80.       MSTK(TOP) = M2   
  81.       NSTK(TOP) = N2   
  82.    13 MN = MSTK(TOP)*NSTK(TOP)            
  83.       CALL WSCAL(MN,SR,SI,STKR(L1),STKI(L1),1)               
  84.       IF (L1.NE.L)     
  85.      $   CALL WCOPY(MN,STKR(L1),STKI(L1),1,STKR(L),STKI(L),1)                   
  86.       GO TO 99         
  87. C                      
  88. C     RIGHT DIVISION   
  89.    20 IF (M2*N2 .EQ. 1) GO TO 21          
  90.       IF (M2 .EQ. N2) FUN = 1             
  91.       IF (M2 .NE. N2) FUN = 4             
  92.       FIN = -1         
  93.       RHS = 2          
  94.       GO TO 99         
  95.    21 SR = STKR(L2)    
  96.       SI = STKI(L2)    
  97.       MN = M*N         
  98.       DO 22 I = 1, MN                     
  99.          LL = L+I-1    
  100.          CALL WDIV(STKR(LL),STKI(LL),SR,SI,STKR(LL),STKI(LL))                   
  101.          IF (ERR .GT. 0) RETURN           
  102.    22 CONTINUE         
  103.       GO TO 99         
  104. C                      
  105. C     LEFT DIVISION    
  106.    25 IF (M*N .EQ. 1) GO TO 26            
  107.       IF (M .EQ. N) FUN = 1               
  108.       IF (M .NE. N) FUN = 4               
  109.       FIN = -2         
  110.       RHS = 2          
  111.       GO TO 99         
  112.    26 SR = STKR(L)     
  113.       SI = STKI(L)     
  114.       MSTK(TOP) = M2   
  115.       NSTK(TOP) = N2   
  116.       MN = M2*N2       
  117.       DO 27 I = 1, MN                     
  118.          LL = L+I-1    
  119.          CALL WDIV(STKR(LL+1),STKI(LL+1),SR,SI,STKR(LL),STKI(LL))               
  120.          IF (ERR .GT. 0) RETURN           
  121.    27 CONTINUE         
  122.       GO TO 99         
  123. C                      
  124. C     POWER            
  125.    30 IF (M2*N2 .NE. 1) CALL ERROR(30)    
  126.       IF (ERR .GT. 0) RETURN              
  127.       IF (M .NE. N) CALL ERROR(20)        
  128.       IF (ERR .GT. 0) RETURN              
  129.       NEXP = IDINT(STKR(L2))              
  130.       IF (STKR(L2) .NE. DFLOAT(NEXP)) GO TO 39               
  131.       IF (STKI(L2) .NE. 0.0D0) GO TO 39   
  132.       IF (NEXP .LT. 2) GO TO 39           
  133.       MN = M*N         
  134.       ERR = L2+MN+N - LSTK(BOT)           
  135.       IF (ERR .GT. 0) CALL ERROR(17)      
  136.       IF (ERR .GT. 0) RETURN              
  137.       CALL WCOPY(MN,STKR(L),STKI(L),1,STKR(L2),STKI(L2),1)   
  138.       L3 = L2+MN       
  139.       DO 36 KEXP = 2, NEXP                
  140.         DO 35 J = 1, N                    
  141.           LS = L+(J-1)*N                  
  142.           CALL WCOPY(N,STKR(LS),STKI(LS),1,STKR(L3),STKI(L3),1)                 
  143.           DO 34 I = 1, N                  
  144.             LS = L2+I-1                   
  145.             LL = L+I-1+(J-1)*N            
  146.             STKR(LL) = WDOTUR(N,STKR(LS),STKI(LS),N,STKR(L3),STKI(L3),1)        
  147.             STKI(LL) = WDOTUI(N,STKR(LS),STKI(LS),N,STKR(L3),STKI(L3),1)        
  148.    34     CONTINUE     
  149.    35   CONTINUE       
  150.    36 CONTINUE         
  151.       GO TO 99         
  152. C                      
  153. C     NONINTEGER OR NONPOSITIVE POWER, USE EIGENVECTORS      
  154.    39 FUN = 2          
  155.       FIN = 0          
  156.       GO TO 99         
  157. C                      
  158. C     ADD OR SUBTRACT SCALAR              
  159.    50 IF (M2 .NE. N2) CALL ERROR(8)       
  160.       IF (ERR .GT. 0) RETURN              
  161.       M = M2           
  162.       N = N2           
  163.       MSTK(TOP) = M    
  164.       NSTK(TOP) = N    
  165.       SR = STKR(L)     
  166.       SI = STKI(L)     
  167.       CALL WCOPY(M*N,STKR(L+1),STKI(L+1),1,STKR(L),STKI(L),1)                   
  168.       GO TO 58         
  169.    52 IF (M .NE. N) CALL ERROR(8)         
  170.       IF (ERR .GT. 0) RETURN              
  171.       SR = STKR(L2)    
  172.       SI = STKI(L2)    
  173.       GO TO 58         
  174.    54 IF (M2 .NE. N2) CALL ERROR(9)       
  175.       IF (ERR .GT. 0) RETURN              
  176.       M = M2           
  177.       N = N2           
  178.       MSTK(TOP) = M    
  179.       NSTK(TOP) = N    
  180.       SR = STKR(L)     
  181.       SI = STKI(L)     
  182.       CALL WCOPY(M*N,STKR(L+1),STKI(L+1),1,STKR(L),STKI(L),1)                   
  183.       CALL WRSCAL(M*N,-1.0D0,STKR(L),STKI(L),1)              
  184.       GO TO 58         
  185.    56 IF (M .NE. N) CALL ERROR(9)         
  186.       IF (ERR .GT. 0) RETURN              
  187.       SR = -STKR(L2)   
  188.       SI = -STKI(L2)   
  189.       GO TO 58         
  190.    58 DO 59 I = 1, N   
  191.          LL = L + (I-1)*(N+1)             
  192.          STKR(LL) = FLOP(STKR(LL)+SR)     
  193.          STKI(LL) = FLOP(STKI(LL)+SI)     
  194.    59 CONTINUE         
  195.       GO TO 99         
  196. C                      
  197. C     COLON            
  198.    60 E2 = STKR(L2)    
  199.       ST = 1.0D0       
  200.       N = 0            
  201.       IF (RHS .LT. 3) GO TO 61            
  202.       ST = STKR(L)     
  203.       TOP = TOP-1      
  204.       L = LSTK(TOP)    
  205.       IF (ST .EQ. 0.0D0) GO TO 63         
  206.    61 E1 = STKR(L)     
  207. C     CHECK FOR CLAUSE                    
  208.       IF (RSTK(PT) .EQ. 3) GO TO 64       
  209.       ERR = L + MAX0(3,IDINT((E2-E1)/ST)) - LSTK(BOT)        
  210.       IF (ERR .GT. 0) CALL ERROR(17)      
  211.       IF (ERR .GT. 0) RETURN              
  212.    62 IF (ST .GT. 0.0D0 .AND. STKR(L) .GT. E2) GO TO 63      
  213.       IF (ST .LT. 0.0D0 .AND. STKR(L) .LT. E2) GO TO 63      
  214.         N = N+1        
  215.         L = L+1        
  216.         STKR(L) = E1 + DFLOAT(N)*ST       
  217.         STKI(L) = 0.0D0                   
  218.         GO TO 62       
  219.    63 NSTK(TOP) = N    
  220.       MSTK(TOP) = 1    
  221.       IF (N .EQ. 0) MSTK(TOP) = 0         
  222.       GO TO 99         
  223. C                      
  224. C     FOR CLAUSE       
  225.    64 STKR(L) = E1     
  226.       STKR(L+1) = ST   
  227.       STKR(L+2) = E2   
  228.       MSTK(TOP) = -3   
  229.       NSTK(TOP) = -1   
  230.       GO TO 99         
  231. C                      
  232. C     ELEMENTWISE OPERATIONS              
  233.    70 OP = OP - DOT    
  234.       IF (M.NE.M2 .OR. N.NE.N2) CALL ERROR(10)               
  235.       IF (ERR .GT. 0) RETURN              
  236.       MN = M*N         
  237.       DO 72 I = 1, MN                     
  238.          J = L+I-1     
  239.          K = L2+I-1    
  240.          IF (OP .EQ. STAR)                
  241.      $      CALL WMUL(STKR(J),STKI(J),STKR(K),STKI(K),STKR(J),STKI(J))          
  242.          IF (OP .EQ. SLASH)               
  243.      $      CALL WDIV(STKR(J),STKI(J),STKR(K),STKI(K),STKR(J),STKI(J))          
  244.          IF (OP .EQ. BSLASH)              
  245.      $      CALL WDIV(STKR(K),STKI(K),STKR(J),STKI(J),STKR(J),STKI(J))          
  246.          IF (ERR .GT. 0) RETURN           
  247.    72 CONTINUE         
  248.       GO TO 99         
  249. C                      
  250. C     KRONECKER        
  251.    80 FIN = OP - 2*DOT - STAR + 11        
  252.       FUN = 6          
  253.       TOP = TOP + 1    
  254.       RHS = 2          
  255.       GO TO 99         
  256. C                      
  257.    99 RETURN           
  258.       END
  259.               
  260.       SUBROUTINE STACKG(ID)               
  261.       INTEGER ID(4)    
  262. C                      
  263. C     GET VARIABLES FROM STORAGE          
  264. C                      
  265.       DOUBLE PRECISION STKR(5005),STKI(5005)                 
  266.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP        
  267.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),HIO,RIO,WIO,RTE,WTE,FE           
  268.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)            
  269.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP          
  270.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,HIO,RIO,WIO,RTE,WTE,FE                  
  271.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN  
  272.       LOGICAL EQID     
  273.       IF (DDT .EQ. 1) WRITE(WTE,100) ID   
  274.   100 FORMAT(1X,'STACKG',4I4)             
  275.       CALL PUTID(IDSTK(1,BOT-1), ID)      
  276.       K = LSIZE+1      
  277.    10 K = K-1          
  278.       IF (.NOT.EQID(IDSTK(1,K), ID)) GO TO 10                
  279.       IF (K .GE. LSIZE-1 .AND. RHS .GT. 0) GO TO 98          
  280.       IF (K .EQ. BOT-1) GO TO 98          
  281.       LK = LSTK(K)     
  282.       IF (RHS .EQ. 1) GO TO 40            
  283.       IF (RHS .EQ. 2) GO TO 60            
  284.       IF (RHS .GT. 2) CALL ERROR(21)      
  285.       IF (ERR .GT. 0) RETURN              
  286.       L = 1            
  287.       IF (TOP .GT. 0) L = LSTK(TOP) + MSTK(TOP)*NSTK(TOP)    
  288.       IF (TOP+1 .GE. BOT) CALL ERROR(18)  
  289.       IF (ERR .GT. 0) RETURN              
  290.       TOP = TOP+1      
  291. C                      
  292. C     LOAD VARIABLE TO TOP OF STACK       
  293.       LSTK(TOP) = L    
  294.       MSTK(TOP) = MSTK(K)                 
  295.       NSTK(TOP) = NSTK(K)                 
  296.       MN = MSTK(K)*NSTK(K)                
  297.       ERR = L+MN - LSTK(BOT)              
  298.       IF (ERR .GT. 0) CALL ERROR(17)      
  299.       IF (ERR .GT. 0) RETURN              
  300. C     IF RAND, MATFN6 GENERATES RANDOM NUMBER                
  301.       IF (K .EQ. LSIZE) GO TO 97          
  302.       CALL WCOPY(MN,STKR(LK),STKI(LK),1,STKR(L),STKI(L),1)   
  303.       GO TO 99         
  304. C                      
  305. C     VECT(ARG)        
  306.    40 IF (MSTK(TOP) .EQ. 0) GO TO 99      
  307.       L = LSTK(TOP)    
  308.       MN = MSTK(TOP)*NSTK(TOP)            
  309.       MNK = MSTK(K)*NSTK(K)               
  310.       IF (MSTK(TOP) .LT. 0) MN = MNK      
  311.       DO 50 I = 1, MN                     
  312.         LL = L+I-1     
  313.         LS = LK+I-1    
  314.         IF (MSTK(TOP) .GT. 0) LS = LK + IDINT(STKR(LL)) - 1  
  315.         IF (LS .LT. LK .OR. LS .GE. LK+MNK) CALL ERROR(21)   
  316.         IF (ERR .GT. 0) RETURN            
  317.         STKR(LL) = STKR(LS)               
  318.         STKI(LL) = STKI(LS)               
  319.    50 CONTINUE         
  320.       MSTK(TOP) = 1    
  321.       NSTK(TOP) = 1    
  322.       IF (MSTK(K) .GT. 1) MSTK(TOP) = MN  
  323.       IF (MSTK(K) .EQ. 1) NSTK(TOP) = MN  
  324.       GO TO 99         
  325. C                      
  326. C     MATRIX(ARG,ARG)                     
  327.    60 TOP = TOP-1      
  328.       L = LSTK(TOP)    
  329.       IF (MSTK(TOP+1) .EQ. 0) MSTK(TOP) = 0                  
  330.       IF (MSTK(TOP) .EQ. 0) GO TO 99      
  331.       L2 = LSTK(TOP+1)                    
  332.       M = MSTK(TOP)*NSTK(TOP)             
  333.       IF (MSTK(TOP) .LT. 0) M = MSTK(K)   
  334.       N = MSTK(TOP+1)*NSTK(TOP+1)         
  335.       IF (MSTK(TOP+1) .LT. 0) N = NSTK(K)                    
  336.       L3 = L2 + N      
  337.       MK = MSTK(K)     
  338.       MNK = MSTK(K)*NSTK(K)               
  339.       DO 70 J = 1, N   
  340.       DO 70 I = 1, M   
  341.         LI = L+I-1     
  342.         IF (MSTK(TOP) .GT. 0) LI = L + IDINT(STKR(LI)) - 1   
  343.         LJ = L2+J-1    
  344.         IF (MSTK(TOP+1) .GT. 0) LJ = L2 + IDINT(STKR(LJ)) - 1                   
  345.         LS = LK + LI-L + (LJ-L2)*MK       
  346.         IF (LS.LT.LK .OR. LS.GE.LK+MNK) CALL ERROR(21)       
  347.         IF (ERR .GT. 0) RETURN            
  348.         LL = L3 + I-1 + (J-1)*M           
  349.         STKR(LL) = STKR(LS)               
  350.         STKI(LL) = STKI(LS)               
  351.    70 CONTINUE         
  352.       MN = M*N         
  353.       CALL WCOPY(MN,STKR(L3),STKI(L3),1,STKR(L),STKI(L),1)   
  354.       MSTK(TOP) = M    
  355.       NSTK(TOP) = N    
  356.       GO TO 99         
  357.    97 FIN = 7          
  358.       FUN = 6          
  359.       RETURN           
  360.    98 FIN = 0          
  361.       RETURN           
  362.    99 FIN = -1         
  363.       FUN = 0          
  364.       RETURN           
  365.       END
  366.               
  367.       SUBROUTINE STACKP(ID)               
  368.       INTEGER ID(4)    
  369. C                      
  370. C     PUT VARIABLES INTO STORAGE          
  371. C                      
  372.       DOUBLE PRECISION STKR(5005),STKI(5005)                 
  373.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP        
  374.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),HIO,RIO,WIO,RTE,WTE,FE           
  375.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)            
  376.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP          
  377.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,HIO,RIO,WIO,RTE,WTE,FE                  
  378.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN  
  379.       LOGICAL EQID     
  380.       INTEGER SEMI     
  381.       DATA SEMI/39/    
  382.       IF (DDT .EQ. 1) WRITE(WTE,100) ID   
  383.   100 FORMAT(1X,'STACKP',4I4)             
  384.       IF (TOP .LE. 0) CALL ERROR(1)       
  385.       IF (ERR .GT. 0) RETURN              
  386.       CALL FUNS(ID)    
  387.       IF (FIN .NE. 0) CALL ERROR(25)      
  388.       IF (ERR .GT. 0) RETURN              
  389.       M = MSTK(TOP)    
  390.       N = NSTK(TOP)    
  391.       IF (M .GT. 0) L = LSTK(TOP)         
  392.       IF (M .LT. 0) CALL ERROR(14)        
  393.       IF (ERR .GT. 0) RETURN              
  394.       IF (M .EQ. 0 .AND. N .NE. 0) GO TO 99                  
  395.       MN = M*N         
  396.       LK = 0           
  397.       MK = 1           
  398.       NK = 0           
  399.       LT = 0           
  400.       MT = 0           
  401.       NT = 0           
  402. C                      
  403. C     DOES VARIABLE ALREADY EXIST         
  404.       CALL PUTID(IDSTK(1,BOT-1),ID)       
  405.       K = LSIZE+1      
  406.    05 K = K-1          
  407.       IF (.NOT.EQID(IDSTK(1,K),ID)) GO TO 05                 
  408.       IF (K .EQ. BOT-1) GO TO 30          
  409.       LK = LSTK(K)     
  410.       MK = MSTK(K)     
  411.       NK = NSTK(K)     
  412.       MNK = MK*NK      
  413.       IF (RHS .EQ. 0) GO TO 20            
  414.       IF (RHS .GT. 2) CALL ERROR(15)      
  415.       IF (ERR .GT. 0) RETURN              
  416.       MT = MK          
  417.       NT = NK          
  418.       LT = L + MN      
  419.       ERR = LT + MNK - LSTK(BOT)          
  420.       IF (ERR .GT. 0) CALL ERROR(17)      
  421.       IF (ERR .GT. 0) RETURN              
  422.       CALL WCOPY(MNK,STKR(LK),STKI(LK),1,STKR(LT),STKI(LT),1)                   
  423. C                      
  424. C     DOES IT FIT      
  425.    20 IF (RHS.EQ.0 .AND. MN.EQ.MNK) GO TO 40                 
  426.       IF (K .GE. LSIZE-3) CALL ERROR(13)  
  427.       IF (ERR .GT. 0) RETURN              
  428. C                      
  429. C     SHIFT STORAGE    
  430.       IF (K .EQ. BOT) GO TO 25            
  431.       LS = LSTK(BOT)   
  432.       LL = LS + MNK    
  433.       CALL WCOPY(LK-LS,STKR(LS),STKI(LS),-1,STKR(LL),STKI(LL),-1)               
  434.       KM1 = K-1        
  435.       DO 24 IB = BOT, KM1                 
  436.         I = BOT+KM1-IB                    
  437.         CALL PUTID(IDSTK(1,I+1),IDSTK(1,I))                  
  438.         MSTK(I+1) = MSTK(I)               
  439.         NSTK(I+1) = NSTK(I)               
  440.         LSTK(I+1) = LSTK(I)+MNK           
  441.    24 CONTINUE         
  442. C                      
  443. C     DESTROY OLD VARIABLE                
  444.    25 BOT = BOT+1      
  445. C                      
  446. C     CREATE NEW VARIABLE                 
  447.    30 IF (MN .EQ. 0) GO TO 99             
  448.       IF (BOT-2 .LE. TOP) CALL ERROR(18)  
  449.       IF (ERR .GT. 0) RETURN              
  450.       K = BOT-1        
  451.       CALL PUTID(IDSTK(1,K), ID)          
  452.       IF (RHS .EQ. 1) GO TO 50            
  453.       IF (RHS .EQ. 2) GO TO 55            
  454. C                      
  455. C     STORE            
  456.    40 IF (K .LT. LSIZE) LSTK(K) = LSTK(K+1) - MN             
  457.       MSTK(K) = M      
  458.       NSTK(K) = N      
  459.       LK = LSTK(K)     
  460.       CALL WCOPY(MN,STKR(L),STKI(L),-1,STKR(LK),STKI(LK),-1) 
  461.       GO TO 90         
  462. C                      
  463. C     VECT(ARG)        
  464.    50 IF (MSTK(TOP-1) .LT. 0) GO TO 59    
  465.       MN1 = 1          
  466.       MN2 = 1          
  467.       L1 = 0           
  468.       L2 = 0           
  469.       IF (N.NE.1 .OR. NK.NE.1) GO TO 52   
  470.       L1 = LSTK(TOP-1)                    
  471.       M1 = MSTK(TOP-1)                    
  472.       MN1 = M1*NSTK(TOP-1)                
  473.       M2 = -1          
  474.       GO TO 60         
  475.    52 IF (M.NE.1 .OR. MK.NE.1) CALL ERROR(15)                
  476.       IF (ERR .GT. 0) RETURN              
  477.       L2 = LSTK(TOP-1)                    
  478.       M2 = MSTK(TOP-1)                    
  479.       MN2 = M2*NSTK(TOP-1)                
  480.       M1 = -1          
  481.       GO TO 60         
  482. C                      
  483. C     MATRIX(ARG,ARG)                     
  484.    55 IF (MSTK(TOP-1).LT.0 .AND. MSTK(TOP-2).LT.0) GO TO 59  
  485.       L2 = LSTK(TOP-1)                    
  486.       M2 = MSTK(TOP-1)                    
  487.       MN2 = M2*NSTK(TOP-1)                
  488.       IF (M2 .LT. 0) MN2 = N              
  489.       L1 = LSTK(TOP-2)                    
  490.       M1 = MSTK(TOP-2)                    
  491.       MN1 = M1*NSTK(TOP-2)                
  492.       IF (M1 .LT. 0) MN1 = M              
  493.       GO TO 60         
  494. C                      
  495.    59 IF (MN .NE. MNK) CALL ERROR(15)     
  496.       IF (ERR .GT. 0) RETURN              
  497.       LK = LSTK(K)     
  498.       CALL WCOPY(MN,STKR(L),STKI(L),-1,STKR(LK),STKI(LK),-1) 
  499.       GO TO 90         
  500. C                      
  501.    60 IF (MN1.NE.M .OR. MN2.NE.N) CALL ERROR(15)             
  502.       IF (ERR .GT. 0) RETURN              
  503.       LL = 1           
  504.       IF (M1 .LT. 0) GO TO 62             
  505.       DO 61 I = 1, MN1                    
  506.          LS = L1+I-1   
  507.          MK = MAX0(MK,IDINT(STKR(LS)))    
  508.          LL = MIN0(LL,IDINT(STKR(LS)))    
  509.    61 CONTINUE         
  510.    62 MK = MAX0(MK,M)                     
  511.       IF (M2 .LT. 0) GO TO 64             
  512.       DO 63 I = 1, MN2                    
  513.          LS = L2+I-1   
  514.          NK = MAX0(NK,IDINT(STKR(LS)))    
  515.          LL = MIN0(LL,IDINT(STKR(LS)))    
  516.    63 CONTINUE         
  517.    64 NK = MAX0(NK,N)                     
  518.       IF (LL .LT. 1) CALL ERROR(21)       
  519.       IF (ERR .GT. 0) RETURN              
  520.       MNK = MK*NK      
  521.       LK = LSTK(K+1) - MNK                
  522.       ERR = LT + MT*NT - LK               
  523.       IF (ERR .GT. 0) CALL ERROR(17)      
  524.       IF (ERR .GT. 0) RETURN              
  525.       LSTK(K) = LK     
  526.       MSTK(K) = MK     
  527.       NSTK(K) = NK     
  528.       CALL WSET(MNK,0.0D0,0.0D0,STKR(LK),STKI(LK),1)         
  529.       IF (NT .LT. 1) GO TO 67             
  530.       DO 66 J = 1, NT                     
  531.          LS = LT+(J-1)*MT                 
  532.          LL = LK+(J-1)*MK                 
  533.          CALL WCOPY(MT,STKR(LS),STKI(LS),-1,STKR(LL),STKI(LL),-1)               
  534.    66 CONTINUE         
  535.    67 DO 68 J = 1, N   
  536.       DO 68 I = 1, M   
  537.         LI = L1+I-1    
  538.         IF (M1 .GT. 0) LI = L1 + IDINT(STKR(LI)) - 1         
  539.         LJ = L2+J-1    
  540.         IF (M2 .GT. 0) LJ = L2 + IDINT(STKR(LJ)) - 1         
  541.         LL = LK+LI-L1+(LJ-L2)*MK          
  542.         LS = L+I-1+(J-1)*M                
  543.         STKR(LL) = STKR(LS)               
  544.         STKI(LL) = STKI(LS)               
  545.    68 CONTINUE         
  546.       GO TO 90         
  547. C                      
  548. C     PRINT IF DESIRED AND POP STACK      
  549.    90 IF (SYM.NE.SEMI .AND. LCT(3).EQ.0) CALL PRINT(ID,K)    
  550.       IF (SYM.EQ.SEMI .AND. LCT(3).EQ.1) CALL PRINT(ID,K)    
  551.       IF (K .EQ. BOT-1) BOT = BOT-1       
  552.    99 IF (M .NE. 0) TOP = TOP - 1 - RHS   
  553.       IF (M .EQ. 0) TOP = TOP - 1         
  554.       RETURN           
  555.       END
  556.               
  557.       SUBROUTINE TERM                     
  558.       DOUBLE PRECISION STKR(5005),STKI(5005)                 
  559.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP        
  560.       INTEGER IDS(4,32),PSTK(32),RSTK(32),PSIZE,PT,PTZ       
  561.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),HIO,RIO,WIO,RTE,WTE,FE           
  562.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)            
  563.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP          
  564.       COMMON /RECU/ IDS,PSTK,RSTK,PSIZE,PT,PTZ               
  565.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,HIO,RIO,WIO,RTE,WTE,FE                  
  566.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN  
  567.       INTEGER R,OP,BSLASH,STAR,SLASH,DOT  
  568.       DATA BSLASH/45/,STAR/43/,SLASH/44/,DOT/47/             
  569.       IF (DDT .EQ. 1) WRITE(WTE,100) PT,RSTK(PT)             
  570.   100 FORMAT(1X,'TERM  ',2I4)             
  571.       R = RSTK(PT)     
  572.       GO TO (99,99,99,99,99,01,01,05,25,99,99,99,99,99,35,99,99,99,99),R        
  573.    01 PT = PT+1        
  574.       RSTK(PT) = 8     
  575. C     *CALL* FACTOR    
  576.       RETURN           
  577.    05 PT = PT-1        
  578.    10 OP = 0           
  579.       IF (SYM .EQ. DOT) OP = DOT          
  580.       IF (SYM .EQ. DOT) CALL GETSYM       
  581.       IF (SYM.EQ.STAR .OR. SYM.EQ.SLASH .OR. SYM.EQ.BSLASH) GO TO 20            
  582.       RETURN           
  583.    20 OP = OP + SYM    
  584.       CALL GETSYM      
  585.       IF (SYM .EQ. DOT) OP = OP + SYM     
  586.       IF (SYM .EQ. DOT) CALL GETSYM       
  587.       PT = PT+1        
  588.       PSTK(PT) = OP    
  589.       RSTK(PT) = 9     
  590. C     *CALL* FACTOR    
  591.       RETURN           
  592.    25 OP = PSTK(PT)    
  593.       PT = PT-1        
  594.       CALL STACK2(OP)                     
  595.       IF (ERR .GT. 0) RETURN              
  596. C     SOME BINARY OPS DONE IN MATFNS      
  597.       IF (FUN .EQ. 0) GO TO 10            
  598.       PT = PT+1        
  599.       RSTK(PT) = 15    
  600. C     *CALL* MATFN     
  601.       RETURN           
  602.    35 PT = PT-1        
  603.       GO TO 10         
  604.    99 CALL ERROR(22)   
  605.       IF (ERR .GT. 0) RETURN              
  606.       RETURN           
  607.       END
  608.               
  609.       SUBROUTINE USER(A,M,N,S,T)          
  610.       DOUBLE PRECISION A(M,N),S,T         
  611. C                      
  612.       INTEGER A3(9)    
  613.       DATA A3 /-149,537,-27,-50,180,-9,-154,546,-25/         
  614.       IF (A(1,1) .NE. 3.0D0) RETURN       
  615.       DO 10 I = 1, 9   
  616.          A(I,1) = DFLOAT(A3(I))           
  617.    10 CONTINUE         
  618.       M = 3            
  619.       N = 3            
  620.       RETURN           
  621.       END 
  622.               
  623.       SUBROUTINE XCHAR(BUF,K)             
  624.       INTEGER BUF(1),K                    
  625. C                      
  626. C     SYSTEM DEPENDENT ROUTINE TO HANDLE SPECIAL CHARACTERS  
  627. C                      
  628. C                      
  629.       INTEGER BACK,MASK                   
  630.       DATA BACK/Z'20202008'/,MASK/Z'000000FF'/               
  631. C                      
  632.       IF (BUF(1) .EQ. BACK) K = -1        
  633.       L = BUF(1) .AND. MASK               
  634.       IF (K .NE. -1) WRITE(6,10) BUF(1),L                    
  635.    10 FORMAT(1X,1H',A1,4H' = ,Z2,' hex is not a MATLAB character.')             
  636.       RETURN           
  637.       END
  638.       SUBROUTINE WGECO(AR,AI,LDA,N,IPVT,RCOND,ZR,ZI)         
  639.       INTEGER LDA,N,IPVT(1)               
  640.       DOUBLE PRECISION AR(LDA,1),AI(LDA,1),ZR(1),ZI(1)       
  641.       DOUBLE PRECISION RCOND              
  642. C                      
  643. C     WGECO FACTORS A DOUBLE-COMPLEX MATRIX BY GAUSSIAN ELIMINATION             
  644. C     AND ESTIMATES THE CONDITION OF THE MATRIX.             
  645. C                      
  646. C     IF  RCOND  IS NOT NEEDED, WGEFA IS SLIGHTLY FASTER.    
  647. C     TO SOLVE  A*X = B , FOLLOW WGECO BY WGESL.             
  648. C     TO COMPUTE  INVERSE(A)*C , FOLLOW WGECO BY WGESL.      
  649. C     TO COMPUTE  DETERMINANT(A) , FOLLOW WGECO BY WGEDI.    
  650. C     TO COMPUTE  INVERSE(A) , FOLLOW WGECO BY WGEDI.        
  651. C                      
  652. C     ON ENTRY         
  653. C                      
  654. C        A       DOUBLE-COMPLEX(LDA, N)   
  655. C                THE MATRIX TO BE FACTORED.                  
  656. C                      
  657. C        LDA     INTEGER                  
  658. C                THE LEADING DIMENSION OF THE ARRAY  A .     
  659. C                      
  660. C        N       INTEGER                  
  661. C                THE ORDER OF THE MATRIX  A .                
  662. C                      
  663. C     ON RETURN        
  664. C                      
  665. C        A       AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS                 
  666. C                WHICH WERE USED TO OBTAIN IT.               
  667. C                THE FACTORIZATION CAN BE WRITTEN  A = L*U  WHERE               
  668. C                L  IS A PRODUCT OF PERMUTATION AND UNIT LOWER                  
  669. C                TRIANGULAR MATRICES AND  U  IS UPPER TRIANGULAR.               
  670. C                      
  671. C        IPVT    INTEGER(N)               
  672. C                AN INTEGER VECTOR OF PIVOT INDICES.         
  673. C                      
  674. C        RCOND   DOUBLE PRECISION         
  675. C                AN ESTIMATE OF THE RECIPROCAL CONDITION OF  A .                
  676. C                FOR THE SYSTEM  A*X = B , RELATIVE PERTURBATIONS               
  677. C                IN  A  AND  B  OF SIZE  EPSILON  MAY CAUSE  
  678. C                RELATIVE PERTURBATIONS IN  X  OF SIZE  EPSILON/RCOND .         
  679. C                IF  RCOND  IS SO SMALL THAT THE LOGICAL EXPRESSION             
  680. C        1.0 + RCOND .EQ. 1.0             
  681. C                IS TRUE, THEN  A  MAY BE SINGULAR TO WORKING                   
  682. C                PRECISION.  IN PARTICULAR,  RCOND  IS ZERO  IF                 
  683. C                EXACT SINGULARITY IS DETECTED OR THE ESTIMATE                  
  684. C                UNDERFLOWS.              
  685. C                      
  686. C        Z       DOUBLE-COMPLEX(N)        
  687. C                A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT.          
  688. C                IF  A  IS CLOSE TO A SINGULAR MATRIX, THEN  Z  IS              
  689. C                AN APPROXIMATE NULL VECTOR IN THE SENSE THAT                   
  690. C                NORM(A*Z) = RCOND*NORM(A)*NORM(Z) .         
  691. C                      
  692. C     LINPACK. THIS VERSION DATED 07/01/79 .                 
  693. C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.              
  694. C                      
  695. C     SUBROUTINES AND FUNCTIONS           
  696. C                      
  697. C     LINPACK WGEFA    
  698. C     BLAS WAXPY,WDOTC,WASUM              
  699. C     FORTRAN DABS,DMAX1                  
  700. C                      
  701. C     INTERNAL VARIABLES                  
  702. C                      
  703.       DOUBLE PRECISION WDOTCR,WDOTCI,EKR,EKI,TR,TI,WKR,WKI,WKMR,WKMI            
  704.       DOUBLE PRECISION ANORM,S,WASUM,SM,YNORM,FLOP           
  705.       INTEGER INFO,J,K,KB,KP1,L           
  706. C                      
  707.       DOUBLE PRECISION ZDUMR,ZDUMI        
  708.       DOUBLE PRECISION CABS1              
  709.       CABS1(ZDUMR,ZDUMI) = DABS(ZDUMR) + DABS(ZDUMI)         
  710. C                      
  711. C     COMPUTE 1-NORM OF A                 
  712. C                      
  713.       ANORM = 0.0D0    
  714.       DO 10 J = 1, N   
  715.          ANORM = DMAX1(ANORM,WASUM(N,AR(1,J),AI(1,J),1))     
  716.    10 CONTINUE         
  717. C                      
  718. C     FACTOR           
  719. C                      
  720.       CALL WGEFA(AR,AI,LDA,N,IPVT,INFO)   
  721. C                      
  722. C     RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) .   
  723. C     ESTIMATE = NORM(Z)/NORM(Y) WHERE  A*Z = Y  AND  CTRANS(A)*Y = E .         
  724. C     CTRANS(A)  IS THE CONJUGATE TRANSPOSE OF A .           
  725. C     THE COMPONENTS OF  E  ARE CHOSEN TO CAUSE MAXIMUM LOCAL                   
  726. C     GROWTH IN THE ELEMENTS OF W  WHERE  CTRANS(U)*W = E .  
  727. C     THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW. 
  728. C                      
  729. C     SOLVE CTRANS(U)*W = E               
  730. C                      
  731.       EKR = 1.0D0      
  732.       EKI = 0.0D0      
  733.       DO 20 J = 1, N   
  734.          ZR(J) = 0.0D0                    
  735.          ZI(J) = 0.0D0                    
  736.    20 CONTINUE         
  737.       DO 110 K = 1, N                     
  738.          CALL WSIGN(EKR,EKI,-ZR(K),-ZI(K),EKR,EKI)           
  739.          IF (CABS1(EKR-ZR(K),EKI-ZI(K))   
  740.      *       .LE. CABS1(AR(K,K),AI(K,K))) GO TO 40           
  741.             S = CABS1(AR(K,K),AI(K,K))    
  742.      *          /CABS1(EKR-ZR(K),EKI-ZI(K))                  
  743.             CALL WRSCAL(N,S,ZR,ZI,1)      
  744.             EKR = S*EKR                   
  745.             EKI = S*EKI                   
  746.    40    CONTINUE      
  747.          WKR = EKR - ZR(K)                
  748.          WKI = EKI - ZI(K)                
  749.          WKMR = -EKR - ZR(K)              
  750.          WKMI = -EKI - ZI(K)              
  751.          S = CABS1(WKR,WKI)               
  752.          SM = CABS1(WKMR,WKMI)            
  753.          IF (CABS1(AR(K,K),AI(K,K)) .EQ. 0.0D0) GO TO 50     
  754.             CALL WDIV(WKR,WKI,AR(K,K),-AI(K,K),WKR,WKI)      
  755.             CALL WDIV(WKMR,WKMI,AR(K,K),-AI(K,K),WKMR,WKMI)  
  756.          GO TO 60      
  757.    50    CONTINUE      
  758.             WKR = 1.0D0                   
  759.             WKI = 0.0D0                   
  760.             WKMR = 1.0D0                  
  761.             WKMI = 0.0D0                  
  762.    60    CONTINUE      
  763.          KP1 = K + 1   
  764.          IF (KP1 .GT. N) GO TO 100        
  765.             DO 70 J = KP1, N              
  766.                CALL WMUL(WKMR,WKMI,AR(K,J),-AI(K,J),TR,TI)   
  767.                SM = FLOP(SM + CABS1(ZR(J)+TR,ZI(J)+TI))      
  768.                CALL WAXPY(1,WKR,WKI,AR(K,J),-AI(K,J),1,      
  769.      $ ZR(J),ZI(J),1)  
  770.                S = FLOP(S + CABS1(ZR(J),ZI(J)))              
  771.    70       CONTINUE   
  772.             IF (S .GE. SM) GO TO 90       
  773.                TR = WKMR - WKR            
  774.                TI = WKMI - WKI            
  775.                WKR = WKMR                 
  776.                WKI = WKMI                 
  777.                DO 80 J = KP1, N           
  778.                   CALL WAXPY(1,TR,TI,AR(K,J),-AI(K,J),1,     
  779.      $    ZR(J),ZI(J),1)                  
  780.    80          CONTINUE                   
  781.    90       CONTINUE   
  782.   100    CONTINUE      
  783.          ZR(K) = WKR   
  784.          ZI(K) = WKI   
  785.   110 CONTINUE         
  786.       S = 1.0D0/WASUM(N,ZR,ZI,1)          
  787.       CALL WRSCAL(N,S,ZR,ZI,1)            
  788. C                      
  789. C     SOLVE CTRANS(L)*Y = W               
  790. C                      
  791.       DO 140 KB = 1, N                    
  792.          K = N + 1 - KB                   
  793.          IF (K .GE. N) GO TO 120          
  794.             ZR(K) = ZR(K)                 
  795.      *            + WDOTCR(N-K,AR(K+1,K),AI(K+1,K),1,ZR(K+1),ZI(K+1),1)         
  796.             ZI(K) = ZI(K)                 
  797.      *            + WDOTCI(N-K,AR(K+1,K),AI(K+1,K),1,ZR(K+1),ZI(K+1),1)         
  798.   120    CONTINUE      
  799.          IF (CABS1(ZR(K),ZI(K)) .LE. 1.0D0) GO TO 130        
  800.             S = 1.0D0/CABS1(ZR(K),ZI(K))  
  801.             CALL WRSCAL(N,S,ZR,ZI,1)      
  802.   130    CONTINUE      
  803.          L = IPVT(K)   
  804.          TR = ZR(L)    
  805.          TI = ZI(L)    
  806.          ZR(L) = ZR(K)                    
  807.          ZI(L) = ZI(K)                    
  808.          ZR(K) = TR    
  809.          ZI(K) = TI    
  810.   140 CONTINUE         
  811.       S = 1.0D0/WASUM(N,ZR,ZI,1)          
  812.       CALL WRSCAL(N,S,ZR,ZI,1)            
  813. C                      
  814.       YNORM = 1.0D0    
  815. C                      
  816. C     SOLVE L*V = Y    
  817. C                      
  818.       DO 160 K = 1, N                     
  819.          L = IPVT(K)   
  820.          TR = ZR(L)    
  821.          TI = ZI(L)    
  822.          ZR(L) = ZR(K)                    
  823.          ZI(L) = ZI(K)                    
  824.          ZR(K) = TR    
  825.          ZI(K) = TI    
  826.          IF (K .LT. N)                    
  827.      *      CALL WAXPY(N-K,TR,TI,AR(K+1,K),AI(K+1,K),1,ZR(K+1),ZI(K+1),         
  828.      *                 1)                 
  829.          IF (CABS1(ZR(K),ZI(K)) .LE. 1.0D0) GO TO 150        
  830.             S = 1.0D0/CABS1(ZR(K),ZI(K))  
  831.             CALL WRSCAL(N,S,ZR,ZI,1)      
  832.             YNORM = S*YNORM               
  833.   150    CONTINUE      
  834.   160 CONTINUE         
  835.       S = 1.0D0/WASUM(N,ZR,ZI,1)          
  836.       CALL WRSCAL(N,S,ZR,ZI,1)            
  837.       YNORM = S*YNORM                     
  838. C                      
  839. C     SOLVE  U*Z = V   
  840. C                      
  841.       DO 200 KB = 1, N                    
  842.          K = N + 1 - KB                   
  843.          IF (CABS1(ZR(K),ZI(K))           
  844.      *       .LE. CABS1(AR(K,K),AI(K,K))) GO TO 170          
  845.             S = CABS1(AR(K,K),AI(K,K))    
  846.      *          /CABS1(ZR(K),ZI(K))       
  847.             CALL WRSCAL(N,S,ZR,ZI,1)      
  848.             YNORM = S*YNORM               
  849.   170    CONTINUE      
  850.          IF (CABS1(AR(K,K),AI(K,K)) .EQ. 0.0D0) GO TO 180    
  851.             CALL WDIV(ZR(K),ZI(K),AR(K,K),AI(K,K),ZR(K),ZI(K))                  
  852.   180    CONTINUE      
  853.          IF (CABS1(AR(K,K),AI(K,K)) .NE. 0.0D0) GO TO 190    
  854.             ZR(K) = 1.0D0                 
  855.             ZI(K) = 0.0D0                 
  856.   190    CONTINUE      
  857.          TR = -ZR(K)   
  858.          TI = -ZI(K)   
  859.          CALL WAXPY(K-1,TR,TI,AR(1,K),AI(1,K),1,ZR(1),ZI(1),1)                  
  860.   200 CONTINUE         
  861. C     MAKE ZNORM = 1.0                    
  862.       S = 1.0D0/WASUM(N,ZR,ZI,1)          
  863.       CALL WRSCAL(N,S,ZR,ZI,1)            
  864.       YNORM = S*YNORM                     
  865. C                      
  866.       IF (ANORM .NE. 0.0D0) RCOND = YNORM/ANORM              
  867.       IF (ANORM .EQ. 0.0D0) RCOND = 0.0D0                    
  868.       RETURN           
  869.       END              
  870.       SUBROUTINE WGEFA(AR,AI,LDA,N,IPVT,INFO)                
  871.       INTEGER LDA,N,IPVT(1),INFO          
  872.       DOUBLE PRECISION AR(LDA,1),AI(LDA,1)                   
  873. C                      
  874. C     WGEFA FACTORS A DOUBLE-COMPLEX MATRIX BY GAUSSIAN ELIMINATION.            
  875. C                      
  876. C     WGEFA IS USUALLY CALLED BY WGECO, BUT IT CAN BE CALLED 
  877. C     DIRECTLY WITH A SAVING IN TIME IF  RCOND  IS NOT NEEDED.                  
  878. C     (TIME FOR WGECO) = (1 + 9/N)*(TIME FOR WGEFA) .        
  879. C                      
  880. C     ON ENTRY         
  881. C                      
  882. C        A       DOUBLE-COMPLEX(LDA, N)   
  883. C                THE MATRIX TO BE FACTORED.                  
  884. C                      
  885. C        LDA     INTEGER                  
  886. C                THE LEADING DIMENSION OF THE ARRAY  A .     
  887. C                      
  888. C        N       INTEGER                  
  889. C                THE ORDER OF THE MATRIX  A .                
  890. C                      
  891. C     ON RETURN        
  892. C                      
  893. C        A       AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS                 
  894. C                WHICH WERE USED TO OBTAIN IT.               
  895. C                THE FACTORIZATION CAN BE WRITTEN  A = L*U  WHERE               
  896. C                L  IS A PRODUCT OF PERMUTATION AND UNIT LOWER                  
  897. C                TRIANGULAR MATRICES AND  U  IS UPPER TRIANGULAR.               
  898. C                      
  899. C        IPVT    INTEGER(N)               
  900. C                AN INTEGER VECTOR OF PIVOT INDICES.         
  901. C                      
  902. C        INFO    INTEGER                  
  903. C                = 0  NORMAL VALUE.       
  904. C                = K  IF  U(K,K) .EQ. 0.0 .  THIS IS NOT AN ERROR               
  905. C  CONDITION FOR THIS SUBROUTINE, BUT IT DOES                
  906. C  INDICATE THAT WGESL OR WGEDI WILL DIVIDE BY ZERO          
  907. C  IF CALLED.  USE  RCOND  IN WGECO FOR A RELIABLE           
  908. C  INDICATION OF SINGULARITY.             
  909. C                      
  910. C     LINPACK. THIS VERSION DATED 07/01/79 .                 
  911. C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.              
  912. C                      
  913. C     SUBROUTINES AND FUNCTIONS           
  914. C                      
  915. C     BLAS WAXPY,WSCAL,IWAMAX             
  916. C     FORTRAN DABS     
  917. C                      
  918. C     INTERNAL VARIABLES                  
  919. C                      
  920.       DOUBLE PRECISION TR,TI              
  921.       INTEGER IWAMAX,J,K,KP1,L,NM1        
  922. C                      
  923.       DOUBLE PRECISION ZDUMR,ZDUMI        
  924.       DOUBLE PRECISION CABS1              
  925.       CABS1(ZDUMR,ZDUMI) = DABS(ZDUMR) + DABS(ZDUMI)         
  926. C                      
  927. C     GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING             
  928. C                      
  929.       INFO = 0         
  930.       NM1 = N - 1      
  931.       IF (NM1 .LT. 1) GO TO 70            
  932.       DO 60 K = 1, NM1                    
  933.          KP1 = K + 1   
  934. C                      
  935. C        FIND L = PIVOT INDEX             
  936. C                      
  937.          L = IWAMAX(N-K+1,AR(K,K),AI(K,K),1) + K - 1         
  938.          IPVT(K) = L   
  939. C                      
  940. C        ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED                  
  941. C                      
  942.          IF (CABS1(AR(L,K),AI(L,K)) .EQ. 0.0D0) GO TO 40     
  943. C                      
  944. C           INTERCHANGE IF NECESSARY      
  945. C                      
  946.             IF (L .EQ. K) GO TO 10        
  947.                TR = AR(L,K)               
  948.                TI = AI(L,K)               
  949.                AR(L,K) = AR(K,K)          
  950.                AI(L,K) = AI(K,K)          
  951.                AR(K,K) = TR               
  952.                AI(K,K) = TI               
  953.    10       CONTINUE   
  954. C                      
  955. C           COMPUTE MULTIPLIERS           
  956. C                      
  957.             CALL WDIV(-1.0D0,0.0D0,AR(K,K),AI(K,K),TR,TI)    
  958.             CALL WSCAL(N-K,TR,TI,AR(K+1,K),AI(K+1,K),1)      
  959. C                      
  960. C           ROW ELIMINATION WITH COLUMN INDEXING             
  961. C                      
  962.             DO 30 J = KP1, N              
  963.                TR = AR(L,J)               
  964.                TI = AI(L,J)               
  965.                IF (L .EQ. K) GO TO 20     
  966.                   AR(L,J) = AR(K,J)       
  967.                   AI(L,J) = AI(K,J)       
  968.                   AR(K,J) = TR            
  969.                   AI(K,J) = TI            
  970.    20          CONTINUE                   
  971.                CALL WAXPY(N-K,TR,TI,AR(K+1,K),AI(K+1,K),1,AR(K+1,J),            
  972.      * AI(K+1,J),1)    
  973.    30       CONTINUE   
  974.          GO TO 50      
  975.    40    CONTINUE      
  976.             INFO = K   
  977.    50    CONTINUE      
  978.    60 CONTINUE         
  979.    70 CONTINUE         
  980.       IPVT(N) = N      
  981.       IF (CABS1(AR(N,N),AI(N,N)) .EQ. 0.0D0) INFO = N        
  982.       RETURN           
  983.       END              
  984.       SUBROUTINE WGESL(AR,AI,LDA,N,IPVT,BR,BI,JOB)           
  985.       INTEGER LDA,N,IPVT(1),JOB           
  986.       DOUBLE PRECISION AR(LDA,1),AI(LDA,1),BR(1),BI(1)       
  987. C                      
  988. C     WGESL SOLVES THE DOUBLE-COMPLEX SYSTEM                 
  989. C     A * X = B  OR  CTRANS(A) * X = B    
  990. C     USING THE FACTORS COMPUTED BY WGECO OR WGEFA.          
  991. C                      
  992. C     ON ENTRY         
  993. C                      
  994. C        A       DOUBLE-COMPLEX(LDA, N)   
  995. C                THE OUTPUT FROM WGECO OR WGEFA.             
  996. C                      
  997. C        LDA     INTEGER                  
  998. C                THE LEADING DIMENSION OF THE ARRAY  A .     
  999. C                      
  1000. C        N       INTEGER                  
  1001. C                THE ORDER OF THE MATRIX  A .                
  1002. C                      
  1003. C        IPVT    INTEGER(N)               
  1004. C                THE PIVOT VECTOR FROM WGECO OR WGEFA.       
  1005. C                      
  1006. C        B       DOUBLE-COMPLEX(N)        
  1007. C                THE RIGHT HAND SIDE VECTOR.                 
  1008. C                      
  1009. C        JOB     INTEGER                  
  1010. C                = 0         TO SOLVE  A*X = B ,             
  1011. C                = NONZERO   TO SOLVE  CTRANS(A)*X = B  WHERE                   
  1012. C         CTRANS(A)  IS THE CONJUGATE TRANSPOSE.             
  1013. C                      
  1014. C     ON RETURN        
  1015. C                      
  1016. C        B       THE SOLUTION VECTOR  X .                    
  1017. C                      
  1018. C     ERROR CONDITION                     
  1019. C                      
  1020. C        A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A           
  1021. C        ZERO ON THE DIAGONAL.  TECHNICALLY THIS INDICATES SINGULARITY          
  1022. C        BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER               
  1023. C        SETTING OF LDA .  IT WILL NOT OCCUR IF THE SUBROUTINES ARE             
  1024. C        CALLED CORRECTLY AND IF WGECO HAS SET RCOND .GT. 0.0                   
  1025. C        OR WGEFA HAS SET INFO .EQ. 0 .   
  1026. C                      
  1027. C     TO COMPUTE  INVERSE(A) * C  WHERE  C  IS A MATRIX      
  1028. C     WITH  P  COLUMNS                    
  1029. C           CALL WGECO(A,LDA,N,IPVT,RCOND,Z)                 
  1030. C           IF (RCOND IS TOO SMALL) GO TO ...                
  1031. C           DO 10 J = 1, P                
  1032. C              CALL WGESL(A,LDA,N,IPVT,C(1,J),0)             
  1033. C        10 CONTINUE   
  1034. C                      
  1035. C     LINPACK. THIS VERSION DATED 07/01/79 .                 
  1036. C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.              
  1037. C                      
  1038. C     SUBROUTINES AND FUNCTIONS           
  1039. C                      
  1040. C     BLAS WAXPY,WDOTC                    
  1041. C                      
  1042. C     INTERNAL VARIABLES                  
  1043. C                      
  1044.       DOUBLE PRECISION WDOTCR,WDOTCI,TR,TI                   
  1045.       INTEGER K,KB,L,NM1                  
  1046. C                      
  1047.       NM1 = N - 1      
  1048.       IF (JOB .NE. 0) GO TO 50            
  1049. C                      
  1050. C        JOB = 0 , SOLVE  A * X = B       
  1051. C        FIRST SOLVE  L*Y = B             
  1052. C                      
  1053.          IF (NM1 .LT. 1) GO TO 30         
  1054.          DO 20 K = 1, NM1                 
  1055.             L = IPVT(K)                   
  1056.             TR = BR(L)                    
  1057.             TI = BI(L)                    
  1058.             IF (L .EQ. K) GO TO 10        
  1059.                BR(L) = BR(K)              
  1060.                BI(L) = BI(K)              
  1061.                BR(K) = TR                 
  1062.                BI(K) = TI                 
  1063.    10       CONTINUE   
  1064.             CALL WAXPY(N-K,TR,TI,AR(K+1,K),AI(K+1,K),1,BR(K+1),BI(K+1),         
  1065.      *                 1)                 
  1066.    20    CONTINUE      
  1067.    30    CONTINUE      
  1068. C                      
  1069. C        NOW SOLVE  U*X = Y               
  1070. C                      
  1071.          DO 40 KB = 1, N                  
  1072.             K = N + 1 - KB                
  1073.             CALL WDIV(BR(K),BI(K),AR(K,K),AI(K,K),BR(K),BI(K))                  
  1074.             TR = -BR(K)                   
  1075.             TI = -BI(K)                   
  1076.             CALL WAXPY(K-1,TR,TI,AR(1,K),AI(1,K),1,BR(1),BI(1),1)               
  1077.    40    CONTINUE      
  1078.       GO TO 100        
  1079.    50 CONTINUE         
  1080. C                      
  1081. C        JOB = NONZERO, SOLVE  CTRANS(A) * X = B             
  1082. C        FIRST SOLVE  CTRANS(U)*Y = B     
  1083. C                      
  1084.          DO 60 K = 1, N                   
  1085.             TR = BR(K) - WDOTCR(K-1,AR(1,K),AI(1,K),1,BR(1),BI(1),1)            
  1086.             TI = BI(K) - WDOTCI(K-1,AR(1,K),AI(1,K),1,BR(1),BI(1),1)            
  1087.             CALL WDIV(TR,TI,AR(K,K),-AI(K,K),BR(K),BI(K))    
  1088.    60    CONTINUE      
  1089. C                      
  1090. C        NOW SOLVE CTRANS(L)*X = Y        
  1091. C                      
  1092.          IF (NM1 .LT. 1) GO TO 90         
  1093.          DO 80 KB = 1, NM1                
  1094.             K = N - KB                    
  1095.             BR(K) = BR(K)                 
  1096.      *            + WDOTCR(N-K,AR(K+1,K),AI(K+1,K),1,BR(K+1),BI(K+1),1)         
  1097.             BI(K) = BI(K)                 
  1098.      *            + WDOTCI(N-K,AR(K+1,K),AI(K+1,K),1,BR(K+1),BI(K+1),1)         
  1099.             L = IPVT(K)                   
  1100.             IF (L .EQ. K) GO TO 70        
  1101.                TR = BR(L)                 
  1102.                TI = BI(L)                 
  1103.                BR(L) = BR(K)              
  1104.                BI(L) = BI(K)              
  1105.                BR(K) = TR                 
  1106.                BI(K) = TI                 
  1107.    70       CONTINUE   
  1108.    80    CONTINUE      
  1109.    90    CONTINUE      
  1110.   100 CONTINUE         
  1111.       RETURN           
  1112.       END              
  1113.       SUBROUTINE WGEDI(AR,AI,LDA,N,IPVT,DETR,DETI,WORKR,WORKI,JOB)              
  1114.       INTEGER LDA,N,IPVT(1),JOB           
  1115.       DOUBLE PRECISION AR(LDA,1),AI(LDA,1),DETR(2),DETI(2),WORKR(1),            
  1116.      *                 WORKI(1)           
  1117. C                      
  1118. C     WGEDI COMPUTES THE DETERMINANT AND INVERSE OF A MATRIX 
  1119. C     USING THE FACTORS COMPUTED BY WGECO OR WGEFA.          
  1120. C                      
  1121. C     ON ENTRY         
  1122. C                      
  1123. C        A       DOUBLE-COMPLEX(LDA, N)   
  1124. C                THE OUTPUT FROM WGECO OR WGEFA.             
  1125. C                      
  1126. C        LDA     INTEGER                  
  1127. C                THE LEADING DIMENSION OF THE ARRAY  A .     
  1128. C                      
  1129. C        N       INTEGER                  
  1130. C                THE ORDER OF THE MATRIX  A .                
  1131. C                      
  1132. C        IPVT    INTEGER(N)               
  1133. C                THE PIVOT VECTOR FROM WGECO OR WGEFA.       
  1134. C                      
  1135. C        WORK    DOUBLE-COMPLEX(N)        
  1136. C                WORK VECTOR.  CONTENTS DESTROYED.           
  1137. C                      
  1138. C        JOB     INTEGER                  
  1139. C                = 11   BOTH DETERMINANT AND INVERSE.        
  1140. C                = 01   INVERSE ONLY.     
  1141. C                = 10   DETERMINANT ONLY.                    
  1142. C                      
  1143. C     ON RETURN        
  1144. C                      
  1145. C        A       INVERSE OF ORIGINAL MATRIX IF REQUESTED.    
  1146. C                OTHERWISE UNCHANGED.     
  1147. C                      
  1148. C        DET     DOUBLE-COMPLEX(2)        
  1149. C                DETERMINANT OF ORIGINAL MATRIX IF REQUESTED.                   
  1150. C                OTHERWISE NOT REFERENCED.                   
  1151. C                DETERMINANT = DET(1) * 10.0**DET(2)         
  1152. C                WITH  1.0 .LE. CABS1(DET(1) .LT. 10.0       
  1153. C                OR  DET(1) .EQ. 0.0 .    
  1154. C                      
  1155. C     ERROR CONDITION                     
  1156. C                      
  1157. C        A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS             
  1158. C        A ZERO ON THE DIAGONAL AND THE INVERSE IS REQUESTED.                   
  1159. C        IT WILL NOT OCCUR IF THE SUBROUTINES ARE CALLED CORRECTLY              
  1160. C        AND IF WGECO HAS SET RCOND .GT. 0.0 OR WGEFA HAS SET                   
  1161. C        INFO .EQ. 0 .                    
  1162. C                      
  1163. C     LINPACK. THIS VERSION DATED 07/01/79 .                 
  1164. C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.              
  1165. C                      
  1166. C     SUBROUTINES AND FUNCTIONS           
  1167. C                      
  1168. C     BLAS WAXPY,WSCAL,WSWAP              
  1169. C     FORTRAN DABS,MOD                    
  1170. C                      
  1171. C     INTERNAL VARIABLES                  
  1172. C                      
  1173.       DOUBLE PRECISION TR,TI              
  1174.       DOUBLE PRECISION TEN                
  1175.       INTEGER I,J,K,KB,KP1,L,NM1          
  1176. C                      
  1177.       DOUBLE PRECISION ZDUMR,ZDUMI        
  1178.       DOUBLE PRECISION CABS1              
  1179.       CABS1(ZDUMR,ZDUMI) = DABS(ZDUMR) + DABS(ZDUMI)         
  1180. C                      
  1181. C     COMPUTE DETERMINANT                 
  1182. C                      
  1183.       IF (JOB/10 .EQ. 0) GO TO 80         
  1184.          DETR(1) = 1.0D0                  
  1185.          DETI(1) = 0.0D0                  
  1186.          DETR(2) = 0.0D0                  
  1187.          DETI(2) = 0.0D0                  
  1188.          TEN = 10.0D0                     
  1189.          DO 60 I = 1, N                   
  1190.             IF (IPVT(I) .EQ. I) GO TO 10  
  1191.                DETR(1) = -DETR(1)         
  1192.                DETI(1) = -DETI(1)         
  1193.    10       CONTINUE   
  1194.             CALL WMUL(AR(I,I),AI(I,I),DETR(1),DETI(1),DETR(1),DETI(1))          
  1195. C           ...EXIT    
  1196. C        ...EXIT       
  1197.             IF (CABS1(DETR(1),DETI(1)) .EQ. 0.0D0) GO TO 70  
  1198.    20       IF (CABS1(DETR(1),DETI(1)) .GE. 1.0D0) GO TO 30  
  1199.                DETR(1) = TEN*DETR(1)      
  1200.                DETI(1) = TEN*DETI(1)      
  1201.                DETR(2) = DETR(2) - 1.0D0  
  1202.                DETI(2) = DETI(2) - 0.0D0  
  1203.             GO TO 20   
  1204.    30       CONTINUE   
  1205.    40       IF (CABS1(DETR(1),DETI(1)) .LT. TEN) GO TO 50    
  1206.                DETR(1) = DETR(1)/TEN      
  1207.                DETI(1) = DETI(1)/TEN      
  1208.                DETR(2) = DETR(2) + 1.0D0  
  1209.                DETI(2) = DETI(2) + 0.0D0  
  1210.             GO TO 40   
  1211.    50       CONTINUE   
  1212.    60    CONTINUE      
  1213.    70    CONTINUE      
  1214.    80 CONTINUE         
  1215. C                      
  1216. C     COMPUTE INVERSE(U)                  
  1217. C                      
  1218.       IF (MOD(JOB,10) .EQ. 0) GO TO 160   
  1219.          DO 110 K = 1, N                  
  1220.             CALL WDIV(1.0D0,0.0D0,AR(K,K),AI(K,K),AR(K,K),AI(K,K))              
  1221.             TR = -AR(K,K)                 
  1222.             TI = -AI(K,K)                 
  1223.             CALL WSCAL(K-1,TR,TI,AR(1,K),AI(1,K),1)          
  1224.             KP1 = K + 1                   
  1225. SHAR_EOF
  1226. #    End of shell archive
  1227. exit 0
  1228. -- 
  1229. Bob Page, U of Lowell CS Dept.  page@swan.ulowell.edu  ulowell!page
  1230. Have five nice days.
  1231.